Title

Subtitle

Author

Affiliation

Ingo Steldermann

Chair of Methods for Model-Based Development in Computational Engineering

Heat transport - explicit time discretization

(supplementary lecture material for the section on the heat equation https://mbd_lectures.pages.rwth-aachen.de/cmm/)

We consider the transport of heat within a vertical material layer, e.g. ice. Let’s start from the heat equation, which reads

\[ \rho c_p \partial_t T(\mathbf{x},t) = \nabla \cdot \left(\kappa \nabla T(\mathbf{x},t)\right). \]

Here, \(T(\mathbf{x},t)\) denotes the temperature field that varies with space \(\mathbf{x}= (x,y,z)\) and time \(t\), \(\rho\) denotes the density, \(c_p\) the specific heat capacity and \(\kappa\) the thermal conductivity.

As stated above, we are interested in the vertical temperature profile, hence restrict ourselves to the \(z\)-dimension. We furthermore divide by \(\rho c_p\) and identify thermal diffusivity as \(\alpha = \kappa/ \rho c_p\), which results in the one-dimensional heat equation:

\[ \partial_t T(z,t) = \alpha \, \partial^2_z T(z,t). \]

Our goal is to solve for the spatio-temporally varying temperature field on the domain

\[[0,z_\mathrm{max}] \times [0,t_\mathrm{end}]\]

subject to initial conditions

\[T(z,0) = T_\mathrm{ini} = \mathrm{const}\]

and boundary conditions

\[T(0,t) = T_\mathrm{surface} \quad \text{and} \quad T(z_\mathrm{max},t) = T_\mathrm{depth}.\]

A simple, straight-forward way to solve this syste is to use finite differences and an explicit time discretization.

Our computational grid is equidistant in both space and time grid \(\{(z_j, t_i)\}\) with

\[z_j \in [0,z_\mathrm{max}], \; 0\leq j \leq M, \quad \text{and} \quad t_i \in [0,t_\mathrm{end}], \; 0\leq i \leq N \]

Resulting grid size is \(\Delta z = z_\mathrm{max}/M\) and time step \(\Delta t\). Ultimately, we are interested in how the vertical temperature profile evolves with time, hence we are interested in

\[T_j^i := T(z_j,t_i)\]

resulting from the initial profile

\[T_j^0 = T_\mathrm{ini} = \mathrm{const}\]

and boundary conditions

\[T_0^i = T_\mathrm{surface} \quad \text{and} \quad T_M^i = T_\mathrm{depth}.\]

Writing the heat equation in finite difference quotients (neglecting higher order terms) yields

\[ \frac{T_j^{i+1}-T_j^i}{\triangle t} = \alpha \frac{T_{j-1}^i-2T_j^i+T_{j+1}^i}{\triangle z^2}, \]

which can be re-organized into an update rule:

\[ T_j^{i+1} = T_j^{i} + \tfrac{\alpha \triangle t}{\triangle z^2} \left(T_{j-1}^{i} - 2 T_{j}^{i} + T_{j+1}^{i}\right) \]

The function euler_expl(…) shows, how this can be implemented.

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
from time import process_time

def euler_expl(alpha, Tini, Tsurf, Tdepth, zmax, M, dt, tend, run_time_plotting = False, version='vectorized'):

    """
    alpha   : thermal diffusivity in [unit]
    Tini    : constant initial condition
    Tsurf   : surface temperature
    Tdepth  : temperature at lower boundary
    zmax    : extent of the computational domain
    M       : number of spatial grid points (M+1)
    dt      : constant time step size
    tend    : end time
    run_time_plotting 
    version : scalar (slow) or vectorized (faster)
    """
 
    # measuring CPU time
    t0 = process_time()  
    
    # determine spatial grid resolution and initialize spatial grid
    dz = zmax/M
    z = np.linspace(0, zmax, M+1)

    # determine number of time steps and initialize time grid
    N = int(round(tend/float(dt)))
    t = np.linspace(0, N*dt, N+1) 

    # make sure dz and dt are compatible with z and t
    dz = z[1] - z[0]
    dt = t[1] - t[0]
    print('Spatial grid resolution: ', dz)
    print('Time step: ', dt)
    
    # determine mesh Fourier number
    F  = alpha * dt / dz**2
    print('Mesh Fourier number: ', round(F,5))
    
    # instantiate temperature vector
    T   = np.zeros(M+1)   
    T_i = np.zeros(M+1)   

    # set initial condition (for all grid points)
    for j in range(0,M+1):
        T_i[j] = Tini

    # check if we want to vizualize during run time
    if run_time_plotting is True:
        fig, ax = plt.subplots()
        plt.gca().invert_yaxis()
        ax.set(xlabel='T(z) [°C]', ylabel='z [m]')
        ax.grid()
        ax.plot(T_i,z, label='temperature profile at t = 0 s')
        fig.canvas.draw()
        plt.pause(1.)

    # loop over all time steps
    for i in range(0, N):
        
        # update all inner points
        if version == 'scalar':
            for j in range(1, M):
                T[j] = T_i[j] + F*(T_i[j-1] - 2*T_i[j] + T_i[j+1])          # <--- here: code was to be completed
        elif version == 'vectorized':
            T[1:M] = T_i[1:M] + F*(T_i[0:M-1] - 2*T_i[1:M] + T_i[2:M+1])    # <--- alternative way: faster
        else:
            raise ValueError('version=%s' % version)

        # insert boundary conditions
        T[0] = Tsurf
        T[M] = Tdepth

        # check if we want to vizualize during run time
        if run_time_plotting is True:
            plt.gca().cla()  # clears current axes
            plt.gca().invert_yaxis()
            ax.set(xlabel='T(z) [°C]', ylabel='z [m]')
            ax.grid()
            ax.plot(T,z, label='temperature profile at t = ' + str(dt*(i+1)) + ' s')
            ax.legend()
            fig.canvas.draw()
            plt.pause(.1)

        # switch variables before next step
        T_i = T
        
        t1 = process_time()
    
    # check if we want to vizualize during run time
    if run_time_plotting is False: # only show final result at tend
        fig, ax = plt.subplots()
        plt.gca().invert_yaxis()
        ax.set(xlabel='T(z) [°C]', ylabel='z [m]')
        ax.grid()
        ax.plot(T,z, label='temperature profile at t = ' + str(dt*N) + ' s')
        ax.legend()
        fig.canvas.draw()
        
        t1 = process_time()
        
    print('Process time [s]', round(t1-t0,3))

    return z,T

Run the algorithm

alpha = 1.203e-6 # thermal diffusivity for ice in m²/s (source: https://de.wikipedia.org/wiki/Temperaturleitfähigkeit)
Tini  = -5.
Tsurf = -10.
Tdepth = 0.
zmax = .1

# numerical parameters
M = 10        # 30 | 10
dt = 20.      # 5 | > 45
tend = 700.   # 750

run_time_plotting = True
version = 'vectorized'

%matplotlib notebook
euler_expl(alpha, Tini, Tsurf, Tdepth, zmax, M, dt, tend, run_time_plotting, version)

Spatial grid resolution:  0.01
Time step:  20.0
Mesh Fourier number:  0.2406
Process time [s] 2.359

(array([0.  , 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1 ]),
 array([-10.        ,  -8.93176928,  -7.88960178,  -6.88960352,
         -5.93177209,  -5.        ,  -4.06822791,  -3.11039648,
         -2.11039822,  -1.06823072,   0.        ]))

Discussion:

The profile relaxes towards the linear profile of the stationary solution. At approx. \(t_\mathrm{end} = [500,750]\,\)s the profile is stationary again. The model is instable if the Fourier number > 0.5; if either \(\triangle t\) or \(\triangle z\) is increased, the other needs to be reduced. For \(M = 10\) the maximum time step \(\triangle t_\mathrm{max,stable} \approx 40\). For increasing \(M\) this value decreases.